An Estimation Method of a Constitutive-law for the Rod Model of 
DNA using Discrete-Structure Simulations 

^Adam R. Hinkle ^Sachin Goyal 

*Harish J. Palanthandalam-Madapusi 

^Theoretical and Applied Mechanics, Cornell University, Ithaca, NY 14853-1503 
* Mechanical and Aerospace Engineering, Syracuse University, Syracuse, NY 13244 

8 June 2009 



Abstract 



The continuum-rod model has emerged as an effi- 
cient tool to describe the long-length-scale structural- 
deformations of DNA which are critical to under- 
standing the nature of many biological processes 
such as gene expression. However, a significant 
challenge in continuum-mechanics-based modeling of 
DNA is to estimate its constitutive law, which follows 
from its interatomic bond-stiffness. Experiments 
and all-atom molecular dynamics (MD) simulations 
have suggested that the constitutive law is non- 
linear and non-homogeneous (sequence-dependent) 
along the length of DNA. In this paper, we present 
an estimation method and a validation study using 
discrete-structure simulations. We consider a simple 
cantilever-rod with an artificially constructed, dis- 
crete lattice-structure which gives rise to a consti- 
tutive law. Large deformations are then simulated. 
An effective constitutive-law is estimated from these 
deformations using inverse methods. Finally, we test 
the estimated law by employing it in the continuum 
rod-model and comparing the simulation results with 
those of discrete-structure simulations under a differ- 
ent cantilever loading-conditions. 



INTRODUCTION 

DNA is a long chain bio-polymer which contains the 
coded, genetic information needed to synthesize pro- 
teins and thus sustain all life [1] . The biological func- 
tions of DNA such as gene expression are known to be 
significantly influenced by its long length-scale struc- 
tural deformations such as looping [2,3], which in 
turn is tied to its chemical make-up — the base-pair 
sequence. For example, the activity of the genes in 
lac-operon in the bacterium E. coli is governed by the 
sequence-dependent looping behavior of its DNA seg- 
ment adjacent to the genes (e.g. refer to [4] and ci- 
tations therein). In fact, this example has become 
a paradigm in understanding looping as a common 
gene-regulatory mechanism. How these long length- 
scale deformations depend on the discrete atomic 
structure molecule is still an open question. 

The DNA molecule itself is found in the nuclei of 
all living cells. It consists of two helices, which wind 
and join together by the hydrogen bonds between 
four bases: adenine (A), cytosine (C), guanine (G), 
and thymine (T). The only bonding which occurs is 
that of A with T (via double hydrogen bonds) and 
C with G (via triple hydrogen bonds) [5,6]. Each 
base in DNA is also attached to a sugar molecule 
and a phosphate molecule. Together, a base, sugar, 
and phosphate are called a nucleotide. Nucleotides 
are arranged in two long strands that form the spi- 
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ral called a double-helix. The structure of the double 
helix is somewhat like a ladder, with the base pairs 
forming the ladder's rungs and the sugar and phos- 
phate molecules forming the vertical sidepieces of the 
ladder as in Figure Q] The sequence of base pairs 
along the helix is the genetic code. 




Figure 1: Three different length scales of the DNA 
molecule. The smallest scale (left) shows double- 
helix structure (sugar-phosphate backbones (orange) 
and base-pairs (blue)). Intermediate scale (middle) 
shows several helical turns of double-stranded DNA 
(ds-DNA) and its helical axis (black). The largest 
scale (right) shows how the helical axis may curve to 
form supercoiling. (Courtesy: Branden and Toozc, 
1999, From Introduction to Protein Structure by Carl 
Branden & John Tooze (reproduced by permission 
of Garland Science/Taylor & Francis, LLC) [7] and 
Lehninger et al. (copied with permission from W.H. 
Freeman) [8]). 

Figure [1] depicts DNA at different length scales. 
The familiar double-helix lies at the smallest scale 
with an approximate diameter of 2 nanometers (nm). 
A complete helical turn extends over nearly 3.6 nm. 
The center images of Figure [T] show the length scale 
of a gene. This solid strand of DNA extends over 
tens to hundreds of helical turns (tens to hundreds of 
nanometers). The base-pair sequence within a gene 
constitutes a chemical code for the production of a 
specific protein via a process known as transcription. 
DNA also needs to make copies of itself via a pro- 
cess known as replication. This is critical when cells 
divide because each new cell needs to have an ex- 
act copy of the DNA present in the old cell. How- 



ever, it is the looping (bending and twisting) at even 
longer length scales which critically influences these 
biological processes of replication and transcription. 
The right portion of Figure Q] illustrates these longer- 
length scale structures which are highly curved and 
twisted topologies called supercoils. The supercoil- 
ing of DNA is also necessary to provide a means for 
these very long molecules to fit within the confines of 
cellular structures, such as the nucleus. 

We are interested in assembling a computational 
framework which can successfully describe and pre- 
dict how DNA undergoes these deformations. Among 
the attempts to model the structural mechanics of 
DNA molecules, few are simple and computationally 
efficient. One such successful approach to modeling 
mechanical behavior of DNA molecules is a contin- 
uum rod model [4,9-12], where elastic properties are 
prescribed and vary with length according to the lo- 
cal base-pair sequence of the molecule. However, a 
key component of the elastic rod model is the con- 
stitutive law (material properties), which is largely 
unknown. There is no clear consensus on its func- 
tional form and how it maps from the base-pair se- 
quence. Although several experiments have provided 
the estimates of linear approximations of average tor- 
sional and bending stiffness by measuring persistence 
lengths [13-16], recent experimental observations and 
model fitting suggest that the constitutive law is non- 
linear [15-19]. The predictions of the continuum rod 
model are extremely sensitive to the constitutive law 
employed and its success hinges upon accurate mod- 
eling of the constitutive law. 

Recently the development of an inverse formula- 
tion of the continuum rod model has been initiated 
to provide an estimation method of the constitutive 
law [17]. In this paper, we extend this method and 
present a validation study using discrete-structure 
simulations of an artificial cantilever-rod. In the 
validation study, we consider a cantilever rod with 
an artificially constructed, discrete, lattice-structure. 
Large deformations in two dimensions are then simu- 
lated. An effective constitutive-law is estimated from 
these deformations using inverse methods. Finally, 
we test the estimated law by employing it in the 
continuum rod and comparing the simulation results 
with those of discrete-structure simulations under a 
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variety of cantilever loading-conditions. Once vali- 
dated, the method can then leverage MD simulations 
of the discrete atomistic structure of DNA and other 
long-chain bio-molecules to estimate their constitu- 
tive laws. 

We begin by first reviewing the mathematical for- 
mulation of the forward rod-model [20] in Section Q] 
and then introduce a general, conceivable form of the 
constitutive law. Then in Section O we describe the 
cantilever example and its validation plan for the esti- 
mation of the constitutive law from some measurable 
outputs. We close by summarizing the conclusions. 




1 The Computational Rod 
Model 

In this section, we review the mathematical formula- 
tion of the forward rod model [20]. The dynamics of 
rod deformation follow from the rigid-body motion 
of its individual cross-sections shown in Figure [2j To 
track the rigid body motion of each cross-section, we 
fix at its mass center a reference frame <Zi(s,i), where 
s is the (unstretched) centerline coordinate and t is 
time. The rigid body motion of the cross-section is 
described by its translational velocity v(s,t) and its 
angular velocity uj(s,t). The gradient of these two 
vector fields along the rod's centerline captures the 
rate of change of the rod's deformation. The curva- 
ture vector, iz(s,t), describes this deformation of the 
rod's centerline (helical axis). The undcformed state 
is denoted by «o(s). Rod deformation in general in- 
volves bending curvature about two axes, shear along 
two axes, torsion, and extension (or compression). 
The position of the cross-section reference frame is 
denoted as R(s,t). The gradient of R(s,t) along s 
is given by r and points along the centerline tan- 
gent t(s,t). In a stress-free state f = t. The stress 
distribution across the cross-section results in a net 
internal (tensile and shear) force /(s, t) and (bending 
and torsional) moment q\s,t). 

The governing equations for the rod dynamics 
are [20]: 
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Figure 2: Dynamical rod model of double-stranded 
DNA on long-length scales. Helical axis of du- 
plex defines the rod centerline which forms a three- 
dimensional space curve located by R(s,t). 
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where m(s) denotes the mass of the rod per unit 
length and the tensor I(s) denotes mass moments 
of inertia per unit length, F(s,t) is the distributed 
external force per unit length and Q(s,t) is the dis- 
tributed external moment per unit length. The par- 
tial derivatives are all relative to the cross-section- 
fixed reference frame dl(s,t). Jl} follows from the 
continuity of R(s,t). |(3J) follows from the continuity 
of the cross-section orientation by constraining the 
curvature (twist) and angular velocity vectors. ^ 
and Q are the Newton-Euler equations for an in- 
finitesimal rod element. 

The stress-strain relationship in the rod integrated 
over its cross section may result in equations of the 
form 

ipi (k, q, f, s) = o. (5) 

([5]) represents the general form of a nonlinear, elas- 
tic constitutive-law. The constitutive law of the rod 
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follows from its interatomic interactions. Therefore 
the varying sequence of base-pairs along the length 
of DNA results in the explicit dependence of ipi on 
s. The result is a non-homogeneous constitutive-law, 
as captured in (0. The basic form of this constitu- 
tive law and its sequence-dependent mapping is an 
open area of research. Current efforts to determine 
the sequence-dependence of the constitutive law al- 
ways begin by assuming a linear approximation to 
© [21,22]. It is important to note, however, that the 
assumption of linearity for the material law of DNA 
has been recently questioned and debated [18,19,23], 
wherein the kink-ability of the molecule indicate non- 
convexity in the constitutive law. 



2 Constitutive Law Estimation 
from a Lattice-Structure 

In this section, we consider a discrete-lattice struc- 
ture of a cantilever in static equilibrium, which is 
fixed at one end and loadings are applied to the free 
end. We construct the discrete-lattice model using 
the commercial, multi-body-dynamics software, Hy- 
perworks. We followed the default system of units in 
Hyperworks, which is [Newton, Millimeter, Kilogram, 
Second] . Figure [3] shows the unit lattice cell where 
point masses are inter-connected by linear springs 
with a spring constant of unity. The point masses 
are arranged to form a unit cube where each mass 
is one unit of distance away from its nearest neigh- 
bors. The springs are connected along each edge of 
the cube and also along the diagonal of each face. 
The cubic cell is repeated thirty times to serve as the 
lumped-parameter model of a cantilever. 

For the static equilibrium of the cantilever, the 
continuum- rod formulation described in Section [T] re- 
duces to a time-independent formulation where all 
partial derivatives with respect to t vanish and v = 
lu = as well. The equations of motion simplify to 
two ordinary vector differential equations in the spa- 
tial domain, s where, 



df 
ds 



+ KX f 



0. 



(7) 



The direction of increasing s is from the clamped 
end to the free end. By the symmetry of the lattice- 
structure, the constitutive law is decoupled in the 
principle directions of bending and torsion and does 
not have any coupling with force f(s,t). Choosing a 
body-fixed frame along the principal directions, ([5]) 
simplifies to three scalar relations, 



ipi (Kt,qi) = 0, 



(8) 



where i — 1,2, or 3 corresponding to the principal 
axes, cii. Specifically, we let a\, a^, and (13 corre- 
spond to the first bending, second bending, and tor- 
sion axes, respectively. 

With this decoupling in place and the assumptions 
of inextensibility and unshearability [20] , we can fur- 
ther simplify our study by looking at deformations 
about one principal axis at a time. Hence we focus 
our further explanation on planar bending. For the 
case of a shear force applied in the ai-direction in- 
plane bending occurs about the a2-axis. The first 
and second components of k and q and the second 
component of / vanish. In this case the governing 
equations of the rod model further reduce to the fol- 
lowing three nonlinear ODEs, 



dq2 
ds 
dfi = 
ds 

dh 
ds 



= -/i, 

■/3«2, 
= /l«2 

and the simplified constitutive law, 



V>2 (^2,92) = 0. 



(9) 
(10) 

(11) 

(12) 



dq 

— + k x q 
ds 



= f x r, 



(6) 



©, (fTU)) , and (fTTj) can be rewritten further in state- 
space form by defining the state vector, x(s), where 
x(s) = [q2{s), fi(s), f 3 {s)] T G R 3 and input u(s) = 
K2- A causal relationship need not exist between in- 
put u and state x. The above state-space equations 
are nonlinear with an unknown algebraic constraint, 
(fT2|) . We note that the state vector x(s) represents 
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the bending moment, shear force, and tension. Sim- 
ilarly, the input u(s) represents the curvature. We 
prescribe the state vector at the free end of the can- 
tilever, i.e. at the final value of s. We convert the con- 
strained final value problem above into a constrained 
initial value problem by substituting s — — s' so that 



dq2 , 

dfi . 
d^ = /3K2 ' 
dh . 

d^ = - flK2 



(13) 
(14) 
(15) 



In terms of the state vector x(s') and input u(s') 
the above equations are of the form 



of this expansion are then found by standard least- 
squares-fitting. 

Once (|T2"j) has been estimated for a certain set of 
loading conditions on the discrete lattice-structure, 
we numerically solve the governing equations given 
by (|16[) for the continuum rod-model, employing 
this estimated constitutive-law. Finally, to vali- 
date the estimated constitutive-law, we apply dif- 
ferent loading conditions to the discrete structure 
and compare these results with corresponding sim- 
ulations of the continuum rod-model using the esti- 
mated constitutive-law. The results of these compar- 
isons determine whether the estimated constitutive- 
law captures the mechanical behavior of the discrete 
structure in its continuum rod-model approximation. 



^ = f(x,u) (16) 

with the constitutive law (|12p as an unknown con- 
straint. 

<j>(x,u)=Q (17) 

The above state-space equations can be solved as 
an initial value problem employing an appropriate 
DAE integrator if the constitutive law is known. We 
use MATLAB Simulink to set up the above equa- 
tions. By contrast, for the estimation of the con- 
stitutive law from some measurable data, we treat 
the inverse method as a state-estimation and input- 
reconstruction problem. For this purpose we employ 
the techniques proposed in [17] which are based on 
filtering and input estimation-algorithms [24-26]. 

In the next section, we illustrate a validation 
scheme for the estimation method in [17]. In particu- 
lar, we assume that the only measurements available 
from the discrete structure simulations are the shape 
of the deformed cantilever and the loads at the free 
end. From the shape, we can compute the curvature 
K 2(s'), which can then be fed as a known input in 
the ODE (fin]) to solve for the state tftj(s / ). To esti- 
mate the functional relationship (| 1 2[) between At2 and 
<72 we express the function, -02, as a sinusoidal-basis- 
function expansion, where the unknown coefficients 




Figure 3: A unit cell of the discrete-lattice model of 
DNA. Eight point masses are displaced a unit dis- 
tance from their nearest-neighbors to form a cube. 
Springs of unit linear stiffness are then attached to 
each mass forming an unstretched configuration. The 
first unit cell, shown here, has the face of the cube 
without springs fixed to resemble the clamp of a can- 
tilever beam. 
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3 Results 



0.1 



We begin by simulating the static equilibrium con- 
figuration of the discrete, cantilever structure with a 
shear force /i = 0.004 as shown in Figure fusing the 
Force Imbalance Method-Type D for the static analy- 
sis in Hyperworks. The system of units was set to the 
default in Hyperworks, which is [Newton, Millimeter, 
Kilogram, Second] . From this configuration data, we 
calculate the curvature K2(s') to use as an input to 
the Simulink model of the continuum rod described 
in the Section [1] The Simulink model calculates the 
bending moment q2{s') according to the continuum- 
rod formulation. Based on these computations, Fig- 
ure [5] shows the estimated constitutive relationship 
between the bending moment and curvature. 

f x =0.004 




Figure 4: Static equilibrium shape of the discrete- 
structure cantilever. Case 1: prescribed shear force, 
but no tension and no bending moment at the free 
end. 

In order to validate the estimated constitutive-law, 
we test it in several different loading environments 
and evaluate if it predicts the same deformation from 
the rod formulation as it does from the discrete- 
structure simulations. Here we present on such val- 
idation test. The loading conditions are shown for 
the validation test in Figure [5] on the deformed equi- 
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Figure 5: Estimated constitutive law from the case 1 
simulation of the discrete structure. 

librium of the discrete structure. In this case we ap- 
ply shear force, compressive force and bending mo- 
ment at the free end. We apply these loading condi- 
tions in the Simulink model of the continuum rod em- 
ploying the estimated constitutive-law shown in Fig- 
ure^ Figure[7]shows the comparison of the centerline 
predicted from discrete-structure simulation (dots) 
and the continuum-rod simulation (solid curve) em- 
ploying the estimated constitutive-law. The closely 
matching shape validates the estimated constitutive- 
law. 

It is likely that the nearly linear result of the esti- 
mated constitutive-law arises from the fact that the 
diagonal springs of the discrete structure are only 
negligibly deformed compared to those parallel to the 
length of the structure. A nonlinear law would be 
expected if deformations of the structure were simu- 
lated such that the diagonal springs experienced large 
deformations. 

The close agreement of the centerline deformations 
of the discrete-lattice structure and the continuum 
rod-model also suggests that a discrete model may 
be an appropriate method to determine a constitu- 
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tive law for DNA. The agreement is noteworthy when 
considering that each simulation is run independently 
with different software and solvers. We also note 
that, while inextensibility and unshearability are as- 
sumed in the forward rod-model, no such assump- 
tions were imposed on the discrete-lattice cantilever. 
The slenderness ratio (length/thickness) of the dis- 
crete structure is 30. Therefore the results suggest 
a confirmation of the assumptions of Kirchhoff Rod 
Theory. 




Figure 6: Static equilibrium shape of the discrete- 
structure cantilever. Case 2: prescribed shear force, 
tension and bending moment at the free end. 



4 Conclusions 

Elastic rod models have recently emerged as efficient 
tools to simulate long length and time scale deforma- 
tions of DNA which are crucial to its biological func- 
tions including gene expression. However, the useful- 
ness of the rod model to simulate DNA is severely 
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Figure 7: Centerline in case 2 loading predicted from 
discrete-structure simulation (dots) and continuum- 
rod simulation (solid curve) employing the estimated 
constitutive-law. 



limited by the lack of knowledge of an accurate con- 
stitutive law for the molecule, which in general, is ex- 
pected to be nonlinear. Discrete structure MD sim- 
ulations of the molecule for small length and time 
scales can provide ample data for the estimation of 
an accurate, constitutive law. However, the past ef- 
forts in mapping the constitutive law from MD simu- 
lations have been limited to linear approximations or 
a few parameter estimations. The challenge has been 
to develop inverse methods to estimate the functional 
form of the constitutive law. 

In this paper, we presented an estimation method 
and a validation study using discrete-structure sim- 
ulations for a simple, but relevant, example moti- 
vated by modeling DNA. We considered a cantilever 
rod with an artificially constructed, discrete, lattice- 
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structure. Large deformations in two dimensions are 
then simulated. An effective constitutive law is esti- 
mated from these deformations using inverse meth- 
ods. Finally, we tested the estimated law by employ- 
ing it in the continuum rod model and comparing 
the simulation results with those of discrete-structure 
simulations under new, cantilever loading-conditions. 
We believe this method can then leverage MD sim- 
ulations of the discrete atomistic structure of DNA 
and other long-chain bio-molcculcs to estimate their 
constitutive laws. 
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